#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <stdbool.h>

#define MAX_ERROR 1e-15
#define MAX_MATRIX_SIZE 100  // 最大矩阵尺寸
#define ROUND_PRECISION 2    // 四舍五入精度（对应Python round(...,2)）

// ========================= 通用工具函数 =========================
// 四舍五入到指定小数位
static double round_to(double val, int precision) {
    double scale = pow(10.0, precision);
    return round(val * scale) / scale;
}

// 数组反转
static void reverse_array(double* arr, int len) {
    for (int i = 0; i < len / 2; i++) {
        double temp = arr[i];
        arr[i] = arr[len - 1 - i];
        arr[len - 1 - i] = temp;
    }
}

// 数组插入元素（头部）
static double* array_insert_head(double* arr, int* len, double val) {
    double* new_arr = (double*)malloc((*len + 1) * sizeof(double));
    new_arr[0] = val;
    memcpy(new_arr + 1, arr, *len * sizeof(double));
    *len += 1;
    free(arr);
    return new_arr;
}

// ========================= 矩阵结构体（模拟Python Matrix类） =========================
typedef struct {
    double** lines;  // 二维数组存储矩阵
    int rows;        // 行数
    int cols;        // 列数
} Matrix;

// 创建矩阵
Matrix* matrix_create(double** lines, int rows, int cols) {
    Matrix* mat = (Matrix*)malloc(sizeof(Matrix));
    mat->rows = rows;
    mat->cols = cols;
    mat->lines = (double**)malloc(rows * sizeof(double*));
    for (int i = 0; i < rows; i++) {
        mat->lines[i] = (double*)malloc(cols * sizeof(double));
        memcpy(mat->lines[i], lines[i], cols * sizeof(double));
    }
    return mat;
}

// 创建空矩阵
Matrix* matrix_create_empty(int rows, int cols) {
    Matrix* mat = (Matrix*)malloc(sizeof(Matrix));
    mat->rows = rows;
    mat->cols = cols;
    mat->lines = (double**)malloc(rows * sizeof(double*));
    for (int i = 0; i < rows; i++) {
        mat->lines[i] = (double*)calloc(cols, sizeof(double));
    }
    return mat;
}

// 释放矩阵内存
void matrix_destroy(Matrix* mat) {
    if (!mat) return;
    for (int i = 0; i < mat->rows; i++) {
        free(mat->lines[i]);
    }
    free(mat->lines);
    free(mat);
}

// 矩阵打印
void matrix_print(const Matrix* mat) {
    for (int i = 0; i < mat->rows; i++) {
        for (int j = 0; j < mat->cols; j++) {
            printf("%.2f ", mat->lines[i][j]);
        }
        printf("\n");
    }
}

// 矩阵拼接（列拼接，对应Python append）
Matrix* matrix_append(const Matrix* mat1, const Matrix* mat2) {
    if (mat1->rows != mat2->rows) {
        fprintf(stderr, "矩阵行数不一致，无法拼接\n");
        return NULL;
    }
    int new_cols = mat1->cols + mat2->cols;
    Matrix* new_mat = matrix_create_empty(mat1->rows, new_cols);
    for (int i = 0; i < mat1->rows; i++) {
        memcpy(new_mat->lines[i], mat1->lines[i], mat1->cols * sizeof(double));
        memcpy(new_mat->lines[i] + mat1->cols, mat2->lines[i], mat2->cols * sizeof(double));
    }
    return new_mat;
}

// 子矩阵提取
Matrix* matrix_sub(const Matrix* mat, int row_start, int row_end, int col_start, int col_end) {
    int rows = row_end - row_start;
    int cols = col_end - col_start;
    Matrix* sub_mat = matrix_create_empty(rows, cols);
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++) {
            sub_mat->lines[i][j] = mat->lines[row_start + i][col_start + j];
        }
    }
    return sub_mat;
}

// 矩阵减法
Matrix* matrix_subtract(const Matrix* mat1, const Matrix* mat2) {
    if (mat1->rows != mat2->rows || mat1->cols != mat2->cols) {
        fprintf(stderr, "矩阵尺寸不一致，无法相减\n");
        return NULL;
    }
    Matrix* res = matrix_create_empty(mat1->rows, mat1->cols);
    for (int i = 0; i < mat1->rows; i++) {
        for (int j = 0; j < mat1->cols; j++) {
            res->lines[i][j] = mat1->lines[i][j] - mat2->lines[i][j];
        }
    }
    return res;
}

// 矩阵乘法
Matrix* matrix_multiply(const Matrix* mat1, const Matrix* mat2) {
    if (mat1->cols != mat2->rows) {
        fprintf(stderr, "矩阵无法相乘（列数≠行数）\n");
        return NULL;
    }
    Matrix* res = matrix_create_empty(mat1->rows, mat2->cols);
    for (int i = 0; i < mat1->rows; i++) {
        for (int j = 0; j < mat2->cols; j++) {
            double sum = 0.0;
            for (int k = 0; k < mat1->cols; k++) {
                sum += mat1->lines[i][k] * mat2->lines[k][j];
            }
            res->lines[i][j] = sum;
        }
    }
    return res;
}

// LUP分解（对应Python lup_decomposition）
int matrix_lup_decomposition(const Matrix* mat, Matrix** L, Matrix** U, Matrix** P) {
    int n = mat->rows;
    if (n != mat->cols) {
        fprintf(stderr, "非方阵无法LUP分解\n");
        return -1;
    }

    // 初始化L/U/P
    *L = matrix_create_empty(n, n);
    *U = matrix_create_empty(n, n);
    *P = matrix_create_empty(n, n);
    double** a = (double**)malloc(n * sizeof(double*));
    int* p = (int*)malloc(n * sizeof(int));

    // 初始化A和置换数组P
    for (int i = 0; i < n; i++) {
        a[i] = (double*)malloc(n * sizeof(double));
        memcpy(a[i], mat->lines[i], n * sizeof(double));
        p[i] = i;
        (*P)->lines[i][i] = 1.0;  // 单位置换矩阵
    }

    for (int i = 0; i < n; i++) {
        // 找主元行（最大绝对值）
        int swap_line = i;
        for (int k = i + 1; k < n; k++) {
            if (fabs(a[k][i]) > fabs(a[swap_line][i])) {
                swap_line = k;
            }
        }

        // 主元为0（无解）
        if (fabs(a[swap_line][i]) <= MAX_ERROR) {
            fprintf(stderr, "矩阵无解（主元为0）\n");
            return -1;
        }

        // 交换行
        if (swap_line != i) {
            double* temp = a[i];
            a[i] = a[swap_line];
            a[swap_line] = temp;

            temp = (*L)->lines[i];
            (*L)->lines[i] = (*L)->lines[swap_line];
            (*L)->lines[swap_line] = temp;

            // 交换置换矩阵P的行
            temp = (*P)->lines[i];
            (*P)->lines[i] = (*P)->lines[swap_line];
            (*P)->lines[swap_line] = temp;

            int p_temp = p[i];
            p[i] = p[swap_line];
            p[swap_line] = p_temp;
        }

        // LU分解核心逻辑
        (*L)->lines[i][i] = 1.0;
        (*U)->lines[i][i] = a[i][i];
        for (int j = i + 1; j < n; j++) {
            (*L)->lines[i][j] = 0.0;
            (*L)->lines[j][i] = round_to(a[j][i] / a[i][i], ROUND_PRECISION);
            (*U)->lines[i][j] = a[i][j];
            (*U)->lines[j][i] = 0.0;

            // 舒尔补计算
            for (int k = i + 1; k < n; k++) {
                a[j][k] = round_to(a[j][k] - (*L)->lines[j][i] * a[i][k], ROUND_PRECISION);
            }
        }
    }

    // 释放临时内存
    for (int i = 0; i < n; i++) {
        free(a[i]);
    }
    free(a);
    free(p);
    return 0;
}

// 解线性方程组（前向替换+后向替换）
double* solve_linear_system(const Matrix* L, const Matrix* U, const double* b, int n, int* out_len) {
    *out_len = n;
    double* y = (double*)malloc(n * sizeof(double));
    double* solution = (double*)malloc(n * sizeof(double));

    // 前向替换解Ly = b
    for (int i = 0; i < n; i++) {
        y[i] = b[i];
        for (int j = i + 1; j < n; j++) {
            y[j] -= L->lines[j][i] * y[i];
        }
    }

    // 后向替换解Ux = y
    for (int i = n - 1; i >= 0; i--) {
        solution[i] = round_to(y[i] / U->lines[i][i], ROUND_PRECISION);
        for (int j = 0; j < i; j++) {
            y[j] -= U->lines[j][i] * solution[i];
        }
    }

    free(y);
    return solution;
}

// 用LUP分解求解线性方程组（对应Python solve_by_lup）
double* matrix_solve_by_lup(const Matrix* mat, const double* values, int val_len, int* out_len) {
    if (mat->rows != val_len) {
        fprintf(stderr, "矩阵行数与值长度不匹配\n");
        return NULL;
    }

    Matrix *L, *U, *P;
    if (matrix_lup_decomposition(mat, &L, &U, &P) != 0) {
        return NULL;
    }

    // 构建b的列矩阵
    double** b_lines = (double**)malloc(val_len * sizeof(double*));
    for (int i = 0; i < val_len; i++) {
        b_lines[i] = (double*)malloc(1 * sizeof(double));
        b_lines[i][0] = values[i];
    }
    Matrix* b_mat = matrix_create(b_lines, val_len, 1);
    for (int i = 0; i < val_len; i++) {
        free(b_lines[i]);
    }
    free(b_lines);

    // 计算Pb
    Matrix* pb_mat = matrix_multiply(P, b_mat);
    double* pb_values = (double*)malloc(val_len * sizeof(double));
    for (int i = 0; i < val_len; i++) {
        pb_values[i] = pb_mat->lines[i][0];
    }

    // 求解
    double* solution = solve_linear_system(L, U, pb_values, val_len, out_len);

    // 释放临时内存
    matrix_destroy(L);
    matrix_destroy(U);
    matrix_destroy(P);
    matrix_destroy(b_mat);
    matrix_destroy(pb_mat);
    free(pb_values);

    return solution;
}

// ========================= 有理函数结构体（模拟Python Rational类） =========================
typedef struct {
    double* numerator;    // 分子多项式系数（反序存储）
    int num_degree;       // 分子次数
    double* denominator;  // 分母多项式系数（反序存储）
    int den_degree;       // 分母次数
} Rational;

// 创建有理函数
Rational* rational_create(double* a, int m, double* b, int n) {
    Rational* r = (Rational*)malloc(sizeof(Rational));
    // 反序存储（对应Python reversed）
    r->num_degree = m;
    r->numerator = (double*)malloc((m + 1) * sizeof(double));
    memcpy(r->numerator, a, (m + 1) * sizeof(double));
    reverse_array(r->numerator, m + 1);

    r->den_degree = n;
    r->denominator = (double*)malloc((n + 1) * sizeof(double));
    memcpy(r->denominator, b, (n + 1) * sizeof(double));
    reverse_array(r->denominator, n + 1);

    return r;
}

// 释放有理函数内存
void rational_destroy(Rational* r) {
    if (!r) return;
    free(r->numerator);
    free(r->denominator);
    free(r);
}

// 多项式计算（Horner法则）
double rational_calc(const double* poly, int len, double val) {
    double result = 0.0;
    for (int i = 0; i < len; i++) {
        result = result * val + poly[i];
    }
    return result;
}

// 多项式转字符串
char* poly_to_str(const double* poly, int len) {
    char* str = (char*)malloc(1024 * sizeof(char));
    str[0] = '\0';
    bool first = true;
    int n = len;

    for (int i = 0; i < len; i++) {
        double coeff = poly[i];
        if (fabs(coeff) < 1e-10) continue;  // 系数为0跳过

        if (!first) {
            strcat(str, "+");
        } else {
            first = false;
        }

        // 系数部分
        if (fabs(coeff) != 1.0 || (n - i - 1) == 0) {
            char coeff_str[64];
            sprintf(coeff_str, "%.6f", coeff);
            strcat(str, coeff_str);
            strcat(str, "*");
        }

        // 次数部分
        int order = n - i - 1;
        if (order > 0) {
            strcat(str, "x");
            if (order > 1) {
                char order_str[16];
                sprintf(order_str, "**%d", order);
                strcat(str, order_str);
            }
        } else {
            strcat(str, "1");
        }
    }

    // 空多项式返回"0"
    if (first) {
        strcpy(str, "0");
    }
    return str;
}

// 有理函数转字符串
char* rational_to_str(const Rational* r) {
    char* num_str = poly_to_str(r->numerator, r->num_degree + 1);
    char* den_str = poly_to_str(r->denominator, r->den_degree + 1);

    char* str = (char*)malloc(2048 * sizeof(char));
    sprintf(str, "(%s)/(%s)", num_str, den_str);

    free(num_str);
    free(den_str);
    return str;
}

// 有理函数求值
double rational_evaluate(const Rational* r, double x) {
    double num_val = rational_calc(r->numerator, r->num_degree + 1, x);
    double den_val = rational_calc(r->denominator, r->den_degree + 1, x);
    return num_val / den_val;
}

// ========================= 帕德逼近核心算法（对应Python pade函数） =========================
Rational* pade_approximation(const double* coefficients, int coeff_len, int m, int n) {
    // 第一步：构建线性方程组求解分母系数
    int mat_rows = n;
    int mat_cols = n;
    double** lines = (double**)malloc(mat_rows * sizeof(double*));

    for (int i = 0; i < n; i++) {
        lines[i] = (double*)malloc(mat_cols * sizeof(double));
        int start = i + 1;
        int end = i + n + 1;
        if (end > coeff_len) {
            fprintf(stderr, "系数长度不足\n");
            return NULL;
        }

        // 提取子数组并反转
        memcpy(lines[i], coefficients + start, n * sizeof(double));
        reverse_array(lines[i], n);
    }

    Matrix* mat = matrix_create(lines, mat_rows, mat_cols);
    for (int i = 0; i < mat_rows; i++) {
        free(lines[i]);
    }
    free(lines);

    // 构建值数组（-coefficients[m+1:m+n+1]）
    int val_len = n;
    double* values = (double*)malloc(val_len * sizeof(double));
    for (int i = 0; i < n; i++) {
        int idx = m + 1 + i;
        if (idx >= coeff_len) {
            fprintf(stderr, "系数长度不足\n");
            return NULL;
        }
        values[i] = -coefficients[idx];
    }

    // 求解分母系数b
    int b_len;
    double* b = matrix_solve_by_lup(mat, values, val_len, &b_len);
    matrix_destroy(mat);
    free(values);

    if (!b) {
        fprintf(stderr, "求解分母系数失败\n");
        return NULL;
    }

    // 插入b0=1
    b = array_insert_head(b, &b_len, 1.0);

    // 第二步：计算分子系数a
    double* a = (double*)malloc((m + 1) * sizeof(double));
    for (int i = 0; i <= m; i++) {
        a[i] = 0.0;
        for (int j = 0; j <= i; j++) {
            if (j >= coeff_len || (i - j) >= b_len) continue;
            a[i] += coefficients[j] * b[i - j];
        }
        a[i] = round_to(a[i], ROUND_PRECISION);
    }

    // 创建有理函数
    Rational* r = rational_create(a, m, b, n - 1);  // b_len = n+1 → 次数n

    // 释放临时内存
    free(a);
    free(b);

    return r;
}

// ========================= 测试代码（sin(x) 3,3阶帕德逼近） =========================
int main() {
    // sin(x)麦克劳林级数系数 [0, 1, 0, -1/6, 0, 1/120, 0]
    double coefficients[] = {0, 1, 0, -1.0/6.0, 0, 1.0/120.0, 0};
    int coeff_len = sizeof(coefficients) / sizeof(double);
    int m = 3, n = 3;

    // 计算帕德逼近
    Rational* r = pade_approximation(coefficients, coeff_len, m, n);
    if (!r) {
        fprintf(stderr, "帕德逼近计算失败\n");
        return -1;
    }

    // 打印结果
    char* r_str = rational_to_str(r);
    printf("帕德逼近结果：%s\n", r_str);

    // 释放内存
    free(r_str);
    rational_destroy(r);

    return 0;
}